Theoretical calculation of the phase behavior of colloidal membranes 
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We formulate a density functional theory that describes the phase behavior of hard rods and depleting poly- 
mers, as realized in recent experiments on suspensions of fd virus and non-adsorbing polymer. The theory 
predicts the relative stability of nematic droplets, stacked smectic columns, and a recently discovered phase 
of isolated monolayers of rods, or colloidal membranes. We find that a minimum rod aspect ratio is required 
for stability of colloidal membranes and that collective protrusion undulations are the dominant effect that sta- 
bilizes this phase. The theoretical predictions are shown to be qualitatively consistent with experimental and 
computational results. 



I. INTRODUCTION 

The relationship between intermolecular interactions and 
structure is fundamental to statistical mechanics and materi- 
als science. Hard particles, which interact solely by steep re- 
pulsive potentials that prohibit overlap, have served as essen- 
tial model systems for understanding this relationship. Stud- 
ies of hard spheres elucidated the structure of liquids |1] and 
3D crystalline phases IJ, [21] and investigations of hard rods 
demonstrated the existence of 3D nematic and smectic phases 
IJ-Q]. Since all accessible hard particle configurations have 
no interparticle interaction energy, these results showed that 
entropy alone can drive the self-assembly of structures with 
long-range order. While the phase diagram of purely repul- 
sive hard rods is well known [60, adding non-adsorbing poly- 
mers introduces attractive interactions between rods through 
the depletion effect, which leads to myriad novel equilibrium 
phases and metastable morphologies that are poorly under- 
stood JT HlOll . Of particular interest from a structural perspec- 
tive, recent experiments on suspensions of monodisperse rod- 
like colloidal viruses and the non-adsorbing polymer Dextran 
demonstrated assembly of 'colloidal membranes' comprised 
of a one rod-length thick monolayer of colloidal rods lUlll . 
This observation has fundamental and practical significance. 
Unlike other examples of entropy-driven assembly which lead 
to long-range order in three dimensions, the colloidal mem- 
branes are self-limited to the thickness of a single rod in one 
dimension and thus are 2D structures. From a practical per- 
spective, equilibrium colloidal membranes could enable man- 
ufacture of inexpensive and easily scalable optoelectronic de- 
vices 1 1211 . 

Previous approaches towards assembly of colloidal mem- 
branes employed chemically heterogeneous rods that mimic 
the amphiphilic nature of lipids which comprise biological 
membranes |13ll . In contrast, the fd molecules involved in as- 
sembly of colloidal membranes are structurally homogeneous, 
suggesting that geometry as well as chemical heterogeneity 
can be used to design molecules that assemble into particu- 
lar structures. Doing so, however, requires a fundamental un- 
derstanding of the forces that conspire to self-limit assembly 
and the relationship between molecular design parameters and 



equilibrium structures. In this article, we therefore construct 
a theoretical description of colloidal rods in the presence of 
depletant molecules. The theory is built on insights from re- 
cent experiments on suspensions of fd virus in Dextran iflHl 
and simulations of hard spherocyUnders in depletant lll4ll . We 
derive a density functional theory for a system of hard cylin- 
ders and depletant. In particular, we use an equation of state 
for hard disks in two dimensions to calculate the equilibrium 
areal density of rods within the membrane, apply free vol- 
ume theory ||15l| to calculate the volume that the membrane 
excludes to polymers, and use a virial expansion to calculate 
rod-rod interactions between nearby membranes. The calcula- 
tions yield predictions for the relationships between osmotic 
pressure, rod aspect ratio, membrane properties, and phase 
behavior. These expressions are derived under the simpli- 
fying assumption that rod orientations are parallel to a fixed 
axis, and are extensively compared to computational and ex- 
perimental results. We find that the theory successfully de- 
scribes the interplay between configurational entropy of rods 
within membranes, depletion interactions, and the resulting 
phase behavior of colloidal rods in depletant. We demonstrate 
that the dominant effect which stabilizes isolated membranes 
arises from the entropic penalty associated with suppression 
of protrusion fluctuations of rods within membranes when the 
membranes stack. While this effect was originally modeled 
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I n ] based on theory describing protrusions of individual rods 
1 1611 . we find that repulsions are driven primarily by collective 
protrusion undulations, consistent with simulations 11411 . 

This paper is organized as follows. We first review the ex- 
periments and simulations on hard rods in the presence of 
depletant molecules in section [III We construct a simplified 
theoretical model for hard rods in depleting polymers and as- 
sociated free energy as a functional of the rod density distri- 
bution in section Hn] We then use the theory to predict the op- 
timal density distributions of rods and hence the system phase 
behavior in section |IV] Under conditions for which mem- 
branes are stable, we analyze the distribution of rods within 
an isolated membrane. We then examine how this distribu- 
tion changes as two membranes approach each other along 
the membrane normals, and thereby evaluate the coupling be- 
tween attractive depletion interactions and repulsive protru- 
sion interactions as a function of membrane separation. Fi- 
nally, predictions for protrusion distributions in isolated and 
interacting membranes, the total interaction potential between 
membranes, and the system phase behavior are compared to 



results of molecular simulations of a rod-depletant suspension 
III4I1 . We find that agreement between the theory and simula- 
tions is quantitative for the nematic/coUoidal membrane co- 
existence osmotic pressure and qualitative for the colloidal 
membrane/smectic-like stacks coexistence pressure. 



II. PREVIOUS EXPERIMENTAL STUDIES OF 
COLLOIDAL MEMBRANES 

We first review recent experiments on suspensions 
of monodisperse rod-like colloidal viruses and the non- 
adsorbing polymer Dextran lUlf (Fig.[T]), which motivate our 
theory, fd viruses alone approximate the behavior of ho- 
mogenous rods interacting with repulsive hard-core interac- 
tions |17|l. The polymer induces an entropy-driven attrac- 
tive (depletion) potential between the rods, the strength and 
range of which can be tuned by changing the polymer con- 
centration and radius of gyration respectively (FiglT]\) H- 
At high polymer concentrations (attraction strength) dilute 
viruses condense into smectic-like stacks of 2D membranes 
(FiglTjD) lUSll . Below a threshold polymer concentration, indi- 
vidual 2D monolayers (membranes) within a smectic filament 
unbind, indicating that the membrane-membrane interaction 
switches from attractive to repulsive 1 1 I1I (FiglTf). The mono- 
layer membranes are stable over months or longer and can be 
many millimeters in diameter. As the polymer concentration 
is decreased further past a second threshold, membranes be- 
come unstable to nematic liquid crystalline droplets or tac- 
toids ([T^) Igl] . While the properties of tactoids and configura- 
tions of rods within them have been explored computationally 
(e.g. iM iOl) and theoretically (e.g. iJbll^BIl ), theoretical 
models of 2D colloidal membranes are lacking. 

Although experiments conclusively demonstrated the exis- 
tence of colloidal membranes, the molecular mechanisms that 
control their stability remained unclear. The depletion inter- 
action that drives lateral association of rods also generates an 
attractive interaction between vertically adjacent membranes. 
For isolated membranes to be stable against stacking, there 
must be a repulsive interaction that overwhelms the deple- 
tion interaction. Experiments on colloidal membranes con- 
taining a low volume fraction of fluorescently labeled rods 
revealed significant protrusions of rods from isolated mem- 
branes, the magnitude of which could be tuned by changing 
the concentration of non-adsorbing polymer In comparison, 
these protrusion fluctuations were suppressed in stacked mem- 
branes Hill . Based on that observation, a model was proposed 
in which the entropy penalty associated with suppressing pro- 
trusion fluctuations of individual rods |lq| leads to repulsive 
interactions membrane-membrane interactions under moder- 
ate osmotic pressure. However, other plausible factors could 
also lead to the observation of isolated membranes, includ- 
ing attractive interactions between virus tips and depletant, 
repulsions due to bending (Helfrich) modes, or kinetic trap- 
ping of metastable membrane intermediates. In Ref. 11411 
we used computer simulations (described in appendix |A] to 
show that collective protrusion effects alone are sufficient to 
produce qualitative agreement with the experimental observa- 
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FIG. 1 . Schematic illustrations and optical micrographs of the self- 
assembled structures observed in suspensions of the filamentous 
virus /rf and non-adsorbing polymer lllll . A) Non-adsorbing poly- 
mer induces effective attractive interactions between rods. B) DIG 
micrograph and schematic of a nematic tactoid formed at low deple- 
tant concentration. C) At intermediate depletant concentrations, rod- 
like viruses condense into macroscopic one rod-length 2D fluid-like 
membranes. D) At high depeltant concentration, membranes stack 
on top of one another, forming smectic filaments. All scale bars are 
5/im. 



tions. The simulations predicted that there is a minimum as- 
pect ratio below which colloidal membranes are never stable; 
this prediction was confirmed by further experiments Ill4ll . 
Here we develop a theory with which we explore the origins 
of the stabilization of colloidal membranes in greater detail. 



III. THEORETICAL MODEL 

In this section we use insights obtained from our previous 
simulations of colloidal membranes 1 14] to derive a tractable 
theoretical model for colloidal membranes. Because the sim- 
ulations suggest that large rod aspect ratios are essential for 
the existence of stable membranes, we cannot model col- 
loidal membranes by directly applying previous theoretical 
approaches used for studying bulk phases of rod-polymer or 
rod-sphere mixtures ll26l430ll . 

Following the simulation model (3, we consider hard rods 
and polymers that can freely interpenetrate other polymers but 
act as hard particles when interacting with rods 131. For sim- 
plicity, rods are represented as cylinders with diameter a and 
length L, and polymers are represented as short cylinders with 
diameter 6 and height ft, = |5. The parameter h is defined 
such that a polymer cylinder has the same volume as a sphere 



of diameter S (the traditional theoretical representation of a 
depletent). Cylinders can interpenetrate one another but ex- 
perience hard-core interactions with rods. Note that because 
we consider an ideal osmolyte (ghost cylinders) we do not 
observe alternating layers of rods and depletents, which were 
described for a model with hard-sphere depletents Bill . 

We focus on conditions relevant to the experiments, where 
rods have large aspect ratios and are immiscible with polymers 
lUin . We showed previously il4ll that under these conditions 
membrane bending modes, which involve deviations of rod 
orientations from the membrane normal [32J], are high-energy 
in comparison to protrusions of rods from the membrane sur- 
face ||33!l on length scales that control membrane-membrane 
stacking. Bending modes can thus can be neglected when 
evaluating phase behavior, and we simplify our calculation by 
constraining rod orientations to be parallel to a fixed direction 
(the z axis). 

We use the density of rods with center of mass at position r, 
p{r), to describe system configurations. To investigate macro- 
scopic membranes we consider periodic boundary conditions 
in the a;y-plane and assume that the density depends only on 
z, p{r) — p{z). The latter simplification follows from con- 
straining rod orientations parallel to z. Then a peak in p{z) 
corresponds to a membrane which is macroscopic in two di- 
mensions (e.g. Fig. |3]below). The width of the peak reflects 
the size of the membrane in the z-direction and thus the extent 
of the protrusion distribution. 

Free volume theory for the free energy. We describe rod- 
rod interactions with a third order virial expansion and rod- 
polymer interactions with the free volume approach presented 
in Ref. II15I1 . adapted to describe the 2-D cross-sections of 
membranes. The results are close to those of a complete third 
order virial expansion, which is lengthy and is not presented 
here. 

In the free volume theory II15I1, for a particular rod den- 
sity distribution p{z) the free energy per unit area, /tot, can be 
written as 



5.y/3/tot =Jdlp{l){lnp{l)X^ - 1) 

dld2p(l)p(2)/(l,2) 

dld2d3p{l)p{2)p{3)f{l, 2)/(l, 3)/(2, 3) 
+ S^yPp, j dz{l-a{z)) (1) 



where Sxy = / dxdy is the total area of the membrane, bold 
numbers are the spatial coordinates, and /(I, 2) is the Mayer 
function between rods. In Eq. [T] the first term is the ideal 
gas free energy and the following two terms are respectively 
the second and third order virial terms for rod-rod interac- 
tions. The second virial term represents the pairwise mutual 
excluded volume interaction between a rod and its neighbors, 
and the third virial term accounts for the mutual exclusion 
among three rods. In the numerical minimization of free en- 
ergy of multi-membranes described below, we found that a 
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FIG. 2. Surface density of rods in isolated membranes. Symbols 
are the results of simulations (appendix IaTi. Lines are predictions 
calculated from the equation of state for a 2-dimensional hard disk 
svstem J34ll as described in the text. The rod lengths, from top to 
bottom, are L = 175, 150, 125, 100, 75 and 50. The diameter of 
a polymer sphere in the simulations is 5 = 1.5 and thus a diameter 
h = lis used for the theoretical cylinders. The dashed line indicates 
the freezing density of a hard disk system, p — 0.88 1351 , 13611 . 



free energy expression with only the second order virial term 
often leads to merging of membranes and unphysically high 
rod densities, thus the three body effects in the third order 
virial term are necessary for physical results. This is to be ex- 
pected, since a second order virial expansion is inaccurate for 
parallel rods ||4|l. 

The last term in Eq. [T] is the free energy due to the vol- 
ume that rods exclude to spheres. The variable a{z) describes 
the free area available to polymers at position z, and de- 
pends on rod densities at any center of mass position from 
which rods can overlap. Specifically, the total density of 
rods that could overlap with a cylinder of height /i at z is 
p*{z) = Jz-L/2-h/2 dzp{z). If we assume that the xy distri- 
bution of rods is not perturbed by polymer, then the fraction of 
free area a{z) can be calculated from p*{z), using the result 
of scaled particle theory for two-dimensions ll37ll . 



a(z) = (! — (/)) cxp(— 277^ ^ JV ^ J V ) 



(2) 



with (/) = 7r(T^p*(z)/4, 7 = ^/(l — 4>), and 77 = 6/ a. Finally 
the excluded volume per unit area for the whole membrane is 
J dz{l — a{z)). By using the scaled particle theory, the effects 
of overlapping excluded volumes of protruding rods are con- 
sidered approximately, and hence the rod-rod correlations are 
partially included. As mentioned above, this approximation 
is equivalent to a third order virial expansion of rod-polymer 
interaction. 

Note that the integration over the z coordinate should be 
restricted to a finite region to avoid divergence, since there 
is always an arbitrarily low but finite concentration of rods 
in the solution, and we will do so in the following numerical 
calculations. 

Since the rod density depends only on z, integration over 
the X and y directions can be carried out analytically to give 



the final form for the free energy as 



tion 



/?/tc 



dzp{z){lnp{z) — 1) 



+ -A I dzip{zi){p\zi + L/2) + p\zi - L/2)) 



+ i^B dzip{zi) I 



dz2p{z2)p\z2+L/2) 



+ -Bj dz,piz^)ipHz,+L/2))^ 
+ Pp, fdz{l~a{z)) 



(3) 



with constants A and B given in appendix |C] and the cumu- 
lative density defined as 



/•z+L/2 

p^z) = / dzp{z) 

Jz-L/2 



(4) 



The details of the calculation are provided in appendix [C] 

The equilibrium rod distribution can be acquired by min- 
imizing /tot with respect to p{z). We will see that Eq. [3] 
can qualitatively describe features of the phase behavior and 
membrane-membrane interactions. However, it does not give 
an accurate prediction of the areal density of rods within the 
membrane, p2d- This limitation is to be expected, since areal 
densities are high and membranes are even crystallized at high 
osmotic pressures according to the simulations. To overcome 
this limitation, we independently relate p2d to the osmotic 
pressure p^, the rod length L and the polymer size h using the 
equation of state for 2D hard disks 113411 . modified to account 
for the extent of the rods in the z direction. The adaptation of 
the equation of state to rods within the membrane is described 
in appendix |B] The equation of state relating p2d to a 2D 
pressure p2d is given in Eqs. IBll to lB4l and the 2D pressure is 
given by p2d ~ (i + /^)ps■ The predicted areal densities are 
compared to simulation results in Fig.|2] 

The equation of state value of the areal density is used as a 
constraint when minimizing the free energy. 



= C 



dzp{z) - mp2d{L,h,p^). 



(5) 



Here the quantity m fixes the total number of rods in the sys- 
tem A^tot as A^tot = rnp2dSxy Under conditions for which iso- 
lated or stacked membranes are stable with respect to the ne- 
matic phase, m will correspond to the number of membranes 
in the system. 

Minimizing /tot thus requires 
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5p{z) Sp{z) 



(6) 



with ( a Lagrange multiplier. Substituting Eqs. [3] and |5] into 
Eq.|6]then results in an integral equation for the rod distribu- 



p{z) =exp(— C) exp(/3ps / dz 



- dp{z) ' 



X cxp{~A{p^z + L/2) + p\z - L/2))) 
X exp{—B dzip{zi)p\zi+ L/2)) 

Jz-L 

xexp(-ii3(pt(, + L/2))2) 



(7) 



The detailed derivation of Eq.|7]is given appendix iDl 

Eq.|7]along with Eq.|5]can be solved numerically to obtain 
the equilibrium rod distribution p{z) for a specified value of 



IV. THEORY RESULTS 

In this section we analyze the behavior predicted by Eqs. [7] 
and |5] For low and moderate osmotic pressures ps, we will see 
that stacks of membranes are thermodynamically unstable; ei- 
ther isolated membranes or nematic configurations (with no 
membranes) are thermodynamically stable. Under these con- 
ditions one can set the number of membranes in Eq. |5] to 
TO = 1 without loss of generahty. Il38ll 



A. Nematic/lsolated membrane phase boundary. 

For low Ps, a flat distribution p{z) — constant is the unique 
solution to Eq. |2l indicating that the nematic state is the ther- 
modynamic equilibrium. (Note that we cannot consider the 
isotropic to nematic transition, which would occur at lower 
Ps, because we have assumed that rods are parallels to the 
z axis.) As Ps increases past a threshold value a stable in- 
homogeneous solution also appears, with a peak in p{z) that 
corresponds to the center of a membrane (Fig. (3]). This solu- 
tion corresponds to an isolated colloidal membrane. The free 
energies of the two solutions are compared to determine the 
equilibrium state. As shown in Fig. IH the predicted coex- 
istence curve for the nematic phase and isolated membranes 
shows remarkable agreement with simulation results for the 
rod lengths considered. As the osmotic pressure increases 
across the spinodal, only the inhomogeneous solution (cor- 
responding to an isolated membrane) remains stable. 

Protrusion distribution. We can further investigate the 
ability of the theoretical model to describe colloidal mem- 
branes by comparing theoretical predictions for the distribu- 
tion of rods within membranes to those measured from sim- 
ulations and a mean field estimate lUol . In the latter ap- 
proach the protrusion of a single rod from a membrane ex- 
posed to depletant osmotic pressure ps incurs a free energy 
/rod = PsAd, with A the cross-sectional area of the rod and 
z the protrusion distance. For uncorrected protrusion sites 
the distribution of protrusions obeys an exponential distribu- 
tion ppi-ot(z) ~ cxp(— PsAz/ZcbT), withppi-ot(z) the density of 
rods with ends located a distance z above the mean surface of 
the membrane. Examples of the numerical solutions for p{z) 
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FIG. 3. Protrusion distribution in a single membrane, z is the dis- 
placement of rods from the center of the membrane. Curves are for 
5 = 1.5, L — 100 and ps = 0.06 (outer lines) or 0.12 (inner lines). 
Solid lines are the simulation results. The dashed lines are the distri- 
butions predicted by EqlT] the solid lines are simulation results, and 
the dotted lines show the scaling expected from an analysis based 
on independent rod protrusions 1101 : p{z) ~ exp(— ps^^), with 
A = 7r(CT + (5)V4. 



are compared to rod distributions measured in simulations of 
isolated membranes and the mean field estimate in Fig|3] We 
see that Eq|2] correctly reproduces the exponential distribu- 
tion at large | z \ (predicted by the mean field estimate) and the 
broadening of the distribution at small \z\. The theory is more 
accurate than the simple mean field estimate 1161] in this con- 
text because the free area accounts for rod-rod correlations 
induced by the overlapping excluded volumes of protruding 
neighbors. 

Although the agreement between theory and simulations is 
good, the predicted distributions are narrower then the simu- 
lation results. We also used the third order virial expansion 
(instead of the free volume approach) to calculate the rod- 
polymer interaction, and obtained similar results. The fact that 
the theory predicts a narrower rod distribution reflects the fact 
that rod-rod correlations are not completely accounted for by 
either the free volume or third order virial expansion calcula- 
tions. The quantitative accuracy likely could be improved by 
going to fourth order in the virial expansion, since graphs in- 
volving polymers (cylinders) that are in the excluded volume 
regions of two neighboring but non-overlapping rods appear 
at this order We also note though that collective protrusions 
introduce long wavelength modes to the membrane, as shown 
by the height-hei ght correlation spectrum (flicker spectrum) 
in Fig. 7 of Ref. Ill4ll . One can potentially account for these 
long wavelength modes using renormalization theory, as in 
Ref. in. 



B. Isolated membrane/smectic phase boundary. 

For osmotic pressures at which membranes are favor- 
able, we determine whether isolated colloidal membranes or 
smectic-like stacks are the thermodynamic minimum by eval- 
uating the rod distribution for the case of two membranes by 
solving Eq. |7] under the constraint of Eq. |5] with ni = 2. 
Numerical solution of Eq. |7] at osmotic pressures above the 
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FIG. 4. (top) Theoretical phase diagram as a function of osmotic 
pressure ps and aspect ratio L. The theoretical phase boundaries, 
calculated as described in the text, are shown as dashed lines. The 
depletant size is 5 = 1.5. (bottom) Phase behavior predicted by sim- 
ulations for the same parameters. Triangles ▲ denote denote parame- 
ters that lead to nematic configurations, -I- symbols correspond to iso- 
lated membranes, and ■ symbols correspond to smectic layers. The 
lower solid line is the theoretical prediction for the nematic/colloidal 
membrane phase boundary, while the upper solid line is fit by eye 
to the the colloidal membrane/smectic phase boundary. The dashed 
line indicates parameter values above which rods crystallized within 
simulated colloidal membranes. Simulation data is from Ref. 11411 . 



nematic-isolated phase boundary yields a stable solution for 
p{z) with two peaks corresponding to two membranes. The 
distance between the peaks depends on the osmotic pressure 
and closely matches simulation results, as shown by the con- 
figuration for Ps = 0.12 and L ~ 100, for which smectic-like 
stacks are thermodynamic ally favorable, in Fig. [5] However, 
there is a finite predicted rod density between the two mem- 
branes, which is likely due to truncating the virial expansion at 
third order, and the predicted membrane widths are narrower 
than those of the simulation, as discussed above. 

Below a certain value of the osmotic pressure, membrane- 
membrane interactions switch from attractive to repulsive, as 
signified by a switch from adjacent to separated peaks in the 
optimal density distribution. However, the exact pressure at 
which this occurs is sensitive to numerical error and depends 
on the region integrated over (which sets the concentration of 
membranes). As noted in Ref. UM the osmotic pressure at 
which smectic-like stacks become stable must depend on the 
concentration of membranes, since the membrane-membrane 
interaction free energy must be sufficiently attractive to com- 
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FIG. 5. Rod density p{z) of two attractive membranes. Tlie the 
daslied line is tlie prediction of Eq.|7]and the solid line is from simu- 
lations. Note that the two peaks are separated by approximately the 
rod length and thus the configuration contains two closely stacked 
membranes. The parameter values are S = 1.5, L = 100, and 
p. = 0.12. 
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FIG. 6. (left) Theoretical free energy of interacting membranes pre- 
dicted by the free energy Eq.[3]with the constraints Eq.|5]and Eq.[8] 
Curves are for parameters 5 = 1.5 and L — 100, with indicated 
values of the osmotic pressures ps. (right) Free energy of interacting 
membranes from simulations using umbrella sampling at the same 
parameters. The solid lines correspond to calculations in which rods 
are constrained parallel to the membrane normals, while the dashed 
line atp.s = 0.06 corresponds to a calculation in which this constraint 
is relaxed. Data is from Ref. u3l. 



pensate the reduction in membrane translational entropy as- 
sociated with stacking. Therefore, to accurately predict the 
colloidal membranes/smectic phase coexistence osmotic pres- 
sure, we next use the theoretical model to calculate the inter- 
action free energy between two membranes. 

Membrane-membrane interactions. To calculate the free 
energy f{d) as a function of the distance between membrane 
centers d we use an approach analogous to umbrella sampling 
II18I1 (this approach is much simpler than performing the pro- 
jection by analytical integration). We augment Eq. |6]with an 
additional constraint on d 



d 
2 



L>o dzzpjz) _ /^^p dzzpjz) 
J dzp{z) J dzp{z) ■ 



(8) 



Following the umbrella sampling procedure 11811 , we imple- 
ment this constraint as a penalty to the free energy: 



/?/t 



penalty 



' L>odzzp{z) d\ 
J.>,dzp{z) 2) 

j( L<odzzp{z) ^ d' 
\Iz<odzp{z) 2 



(9) 



with fc > an adjustable constant. We then numerically mini- 
mize /tot + /penalty to obtain the optimal density distribution 
p{z; d) under the constraint Eq. [8] Finally, the interaction 
free energy f{d) is obtained by subtracting the penalty term: 
.f{d) = .ftot{p{z;d))/mp2d- 

Examples of interaction free energies are shown in Fig. |6] 
(top). If we compare these theoretical predictions to simula- 
tion umbrella sampling results 114!] in Fig.|6](bottom), we see 
that the theoretical f{d) curves qualitatively agree with sim- 
ulation results, but that the theoretical calculations have shal- 
lower attractive basins. Like the discrepancy between the the- 
oretical and computational protrusion distributions, this quan- 
titative difference may be due to the fact that the rod-rod cor- 
relations are not fully accounted for 
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FIG. 7. Theoretical calculation of the total free energy of attractive 
basin F, defined by Eq. [TO]with M = 10"*, for rod lengths L = 
40, 75, 100 and 150 (solid lines, from left to right). The horizontal 
dashed line is Fq, defined by Eq.QTjwith pmVo = 10^*. 



Following Ref. |[l4il, smectic layers are thermodynamically 
favorable at finite membrane concentration p,^ if the total free 
energy of the attractive basin in the membrane-membrane in- 
teraction potential satisfies 



F < Fq = kBTlnp^vo 



(10) 



with 



exp(- 



-I3F) = f 
J f 



/(s)<0 



dsexp(-2/3M/(s)) (11) 



with M the number of rods in one membrane, and uq a stan- 
dard state volume. We roughly estimate M = 10** and 
Pm^'o = 10^* from the experimental conditions; the loca- 
tion of the phase boundary is not sensitive to the value of 
Prn^'o- Fig. [T] shows the theoretical values for the free en- 
ergy of attractive basins F at a number of rod lengths. As 
expected, F becomes more favorable as the osmotic pressure 
Ps increases. The theoretical F curves cross k^TXwp^VQ near 
Ps ~ 0.07, which is close to the simulation results for deple- 
tant size 6 = 1.5. 

As shown in Fig. H] the theoretical isolated-smectic phase 
boundary shows reasonable agreement with the simulated 



phase boundary for L < 100. Furthermore, both methods 
predict a similar threshold value of the aspect ratio, L w 30, 
below which the system transitions directly from nematic con- 
figurations to smectic-like stacks of membranes. This pre- 
diction was confirmed by experiments 1 14il in which the os- 
motic pressure was varied by controlling the concentration of 
non-adsorbing polymer and the depletant size was varied by 
changing polymer radius of gyration. 

Notably, the theory does not reproduce the decrease in the 
transition osmotic pressure at large rod lengths seen in the 
simulations. As described in Ref. 1114(1 . this trend results from 
crystallization of rods within membranes in the simulations. 
Under the simulation conditions, membrane crystallization 
decreases the interaction free energy due to protrusions and 
thus lowers the isolated-smectic coexistence osmotic pressure. 
In contrast, the theoretical interaction free energy f{d) as- 
sumes a disordered distribution of rod positions within the 
plane of the membrane and thus does not allow for membrane 
crystallization. It is worth noting that simulations in which the 
constraint on parallel rod orientations was relaxed required 
larger rod aspect ratios and/or higher osmotic pressures for 
crystallization of rods within membranes. Given that observa- 
tion and the fact that theory and simulation agree for L < 100, 
we chose not to extend the theoretical model to allow for crys- 
tallization. 



cation, it will also be desirable to allow for rod orientational 
fluctuations, to determine their effect on the locations of phase 
boundaries and to determine the effects of membrane bending 
modes on membrane-membrane interactions. As discussed 
earlier, these modes are expected to be of limited importance 
under the experimental conditions due to the high membrane 
bending modulus, but will become more important as the rod 
length is decreased. Finally, it would be useful to include ef- 
fects of semi-flexibility. It was shown through scaling argu- 
ments in Ref. lll4ll that semi-flexibility renormalizes the in- 
teractions between rods in the membrane leading to a smaller 
equilibrium areal density; in addition, semiflexible rods will 
behave as if polydisperse in length. 
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V. CONCLUSIONS AND OUTLOOK 

In summary, we have presented a theoretical model that 
represents hard rods in the presence of depletant molecules 
such as non-adsorbing polymer. The free energy is con- 
structed by using the free volume theory for depletant-rod in- 
teractions, and a third order virial expansion for rod-rod in- 
teractions, with the equation of state for a hard disk system 
to constrain the areal rod density. The predicted nematic- 
membrane phase boundary shows reasonable agreement with 
simulation results, and the predicted isolated-smectic phase 
boundary qualitatively agrees with simulation results for rod 
length L < lOOcr. The predicted phase boundaries establish 
that there is a critical rod length L ^ SOcr, below which iso- 
lated colloidal membrane can not be formed, for the given 
depletant size S = l.Scr. 

The theoretical calculations enable systematic identifica- 
tion of the factors that control system phase behavior. In par- 
ticular, the theoretical results demonstrate that correlations be- 
tween protruding rods significantly enhance protrusion fluctu- 
ations and thereby profoundly affect membrane morphologies 
and the range of interactions between membranes. This effect 
is emphasized by comparison of the theoretically predicted 
protrusion distributions with those of a simpler theory that ne- 
glects correlations between neighboring rods. Evaluation of 
the theoretical predictions at different orders of the virial ex- 
pansion further demonstrate that the terms up to third order 
that we have considered are essential for physical predictions. 
We speculate that including fourth order terms would lead to 
more quantitative agreement with the simulation results below 
the membrane crystallization point. At this level of sophisti- 



Appendix A: Simulations of colloidal membranes 



In this section we describe the computer simulations whose 
results we have compared with the theory predictions. The 
simulations describe the equilibrium phase behavior for a 
model of hard rods and depletant molecules in the absence 
of any attractive interactions between rod ends and depletant 
III4I1 . The rods are represented as hard spherocylinders with 
diameter a and length L. The non-adsorbing polymer (deple- 
tant) is modeled with ghost spheres |4C|| of diameter 6, which 
act as hard spheres when interacting with rods but can freely 
interpenetrate one another Compared with an effective pair 
potential approach ll4ll - l43ll . this model accounts for multi-rod 
interactions induced by polymers. Simulation results are re- 
ported with the following units: a is the unit of length, ksT is 
the unit of energy, and ksTa^'^ is the unit of n-dimensional 
pressure (n = 2 or 3). As noted above, tilting of rods does 
not qualitatively affect membrane-membrane interactions and 
thus most free energy simulations were performed with rod 
orientations constrained parallel to the z axis of the simula- 
tion box. 

Examples of membrane-membrane interaction free ener- 
gies calculated by umbrella sampling II18I1 are shown in Fig. |6] 
(bottom). For p^ = 0.06 the free energy calculated under the 
orientational constraint is compared to the free energy calcu- 
lated with the constraint relaxed. The phase behavior of the 
computational model system was predicted as a function of 
the depletant osmotic pressure p,, the rod aspect ratio L, and 
the size of the depletant ghost spheres 5, which corresponds 
to the polymer radius of gyration. Fig. |4] (bottom) shows a 
cross-section of the phase diagram in terms of ps and L. 
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Appendix B: Areal rod densities 

As noted in the main text, the virial expansion Eq.[T] can- 
not accurately predict the areal densities of rods p2d within a 
membrane because areal densities are high under conditions 
for which membranes are stable, and membranes crystallize 
under high osmotic pressures. 

We therefore independently obtain p2d by noting that, under 
the assumption that rods are parallel to z direction, the cross 
section of a colloidal membrane can be considered as a 2D 
system of hard disks. We ignore the small density of polymer 
that is actually inside the membrane, and the osmotic pressure 
Ps determines the two-dimensional pressure felt by the hard 
disks (rods). The equation of state for hard disk systems has 
been studied extensively y^Bj-fial- Here we use the global 
equation of state given by Luding Il34ll which characterizes 
both the liquid and crystalline phases. The equation of state is 
given by 



PPld/Pld -1 = Pi+ m(l^)(Pden,se - Pi) 



(Bl) 



where v ~ P2d7r 0-^/4 is the areal fraction. P4 is the low den- 
sity result 



with 



9i{^) 



2vgi[v) 



1 - 7i^/16 v^/W 



(1 



[l-uY 



and Pdense IS the high density result 

Co 



-'dense 



hsil^mi^x -v) -1 



(B2) 



(B3) 



(B4) 



with i/max = 7r/(2-\/3) the maximum areal fraction, h^{x) = 
1 + cix + c^x^ a fit polynomial, and constants cp = 1.8137, 
ci = -0.04 and C3 = 3.25. Full details are in Q. Eq.lBT] 
gives the pressure p2d as a function of the density, but the den- 
sity can be numerically inverted for a given pressure to give: 



P2d = P2d(P2d) 



(B5) 



The 2D pressure p2d is the result of polymer osmotic pres- 
sure acting laterally to the membrane. The excluded volume 
per rod in the membrane is approximately Wex ~ (^ + h)/ p2d, 
where we neglect rod protrusions. Mechanical equilibrium 
then requires the 2D pressure to be 



P2d = -Pid^T—Ps W (i + h)ps 



'dp: 



(B6) 



'2d 



The areal density of rods p2d{L, h,Pi) is then acquired from 
Eq. lBSI and Eq. lB6l for a given L, h, and p^. 

Predicted values of p2d are compared to simulation results 
in Fig.|2] The simulation results of y02d are obtained from sta- 
ble membranes of 256 rods; areal densities are not sensitive 
to system size. We find that the difference between simulated 
and predicted areal densities is within 5%, and estimated val- 
ues are always lower than the simulation results. The differ- 
ence can be attributed to the rough estimation of the excluded 



volume per rod and the neglect of rod protrusions, as well as 
the accuracy of Eq. lBni34ll . The sharp increase of /92d, starting 
near p2d ~ 0.88, identifies the transition from liquid phase to 
crystal phase (identified from the two-dimensional radial dis- 
tribution function of rods q(r) 114)1 ). with p2d ~ 0.88 the hard 
disk freezing point OSl 13611 . 



Appendix C: Integration of free energy 

Because rods are required to be parallel to z direction, the 
free energy expression Eq. [T] can be greatly simplified. The 
rod-rod Mayer function /(I, 2) = exp(— /?C/(1, 2)) - 1 can 
be separated into to the Mayer function in z direction and the 
Mayer function in x — y plane 114711 . 



with 



/.(I, 2) 



/(1,2) = -/,(1,2)/,,(1,2) 



-1, |zi - Z2I < L 
0, otherwise 



(CI) 



/.,(i,2) = |-^' i-^--2r + iy^-y2r<'^' (C2) 

•'^^ ' ' 10, otherwise ^ ^ 



Since the rod density depends only on z, the integrations in 
Eq.[T]can be separated as well. We have. 



/?/= dzpiz)iliipiz)-l) 



- l^ J lf[dz,p{zA Ml,2) 

-\b j [\{d^^Pi^i)\ /.(l,2)/,(l,3)/,(2,3) 



+ f3ps I dz{l — a{z)) 



(C3) 



with 



^^-^J(l[ dx^dy, j /,.,(!, 2) 



(C4) 



and 



B =- ^ j \\{dxdyi\ f.y{l,2)Uy{\,^)f^y{2,Z) 



4(2 3v3 



(C5) 
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The integrations over z can be further simplified, 

dzip(zi) / dz2p{z2) 

Jzi-L 

dzipizi){p\zi + L/2) + pt(zi - L/2)) (C6) 



and 



l[dz,p{z,)] M1,2)M1,3)M2,3) 



\i=l 



pZi pZ2+L 

dzip{zi) / dz2p{z2) / dz3p(z3) 

J Z\—L J Z2 

{Z2 ^ 23) 

/•zi+L rzi+L 

dzip{zi) / dz2p{z2) / dzspizs) 

J Z\ J Z\ 

-2 / dzi/9(^l) / dZ2/0(^2)/5t(2;2 + ^/2) 

J Jzi-L 



dz^p{zi){p\zi+L/2)f 



(C7) 



with ^^(z) defined in Eq.|4] The final expression for the free 



energy, Eq. [3] is then acquired by substituting equations from 
Eq.l^to Eq.|C2]into Eq.|C3] 



Appendix D: Free energy minimization 



The two terms in Eq.|6]are calculated as, 

:lnp(z) + A{p\z + L/2) + p\z- L/2)) 



SPftot 



Sp{z) 



b[ dz2p{z2)p\z2 + L/2) 

Jz-L 



+ ^B{pHz,+L/2))' 



PPs / dzi 



da{zi) 
dp{z) 



and 



5C 

J^) 



= 1 



(Dl) 



(D2) 



Note that the factor 1/2 in front of A and the factor 1/3 in 
front of B are canceled because p appears multiple times in 
the corresponding terms (see Eq.[T]). Substituting these results 
into Eq. |6] the resulting integral equation of p is acquired as 
Eq. |7] When p is numerically solved, the constraint Eq. |5]is 
applied through Q in each iteration step. 
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